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Abstract 

Dynamics of spiral waves in perturbed, e.g. slightly inhomogeneous or subject to a small 
periodic external force, two-dimensional autowave media can be described asymptotically in 
terms of Aristotelean dynamics, so that the velocities of the spiral wave drift in space and 
time are proportional to the forces caused by the perturbation. The forces are defined as a 
convolution of the perturbation with the spirals Response Functions, which are eigenfunctions 
of the adjoint linearised problem. In this paper we find numerically the Response Functions 
of a spiral wave solution in the classic excitable FitzHugh-Nagumo model, and show that they 
are effectively localised in the vicinity of the spiral core. 

1 Introduction 

Autowaves are nonlinear waves observed in spatially distributed media of physical, chemical, and 
biological nature, where wave propagation is supported by a source of energy stored in the medium. 
In a two-dimensional autowave medium there may exist autowave vortices appearing as rotating 
spiral waves and thus acting as sources of periodic waves. Their existence is not due to singularities 
in the medium but is determined only by development from initial conditions. In a slightly 
perturbed medium, e.g. spatially inhomogeneous or subject to time-dependent external forcing, a 
spiral wave drifts, i.e. its core location and frequency change with time (Biktashev 2004, Biktasheva 
2000, Biktasheva, Elkin & Biktashev 1999, Fast & Pcrtsov 1990, Pcrtsov & Ermakova 1988). 

While the hypothesis of re-entry of excitation underlying cardiac arrhythmias belongs to the 
beginning of the twentieth century, e.g. (Mines 1913), the first direct experimental observation of 
spiral waves was reported in 1960s in a chemical oscillatory medium, the Belousov-Zhabotinsky 
(BZ) reaction (Zhabotinsky & Zaikin 1971). That triggered a huge amount of interest and activity 
in the area. Soon after that spiral waves were observed in a rabbit ventricular tissue (Allessie, Bonk 
& Schopman 1973), and later in a variety of other spatially distributed active systems: in chick 
retina (Gorelova & Bures 1983), colonies of social amoebae (Alcantara & Monk 1974), cytoplasm of 
single oocytes (Lechleiter, Girard, Peralta & Clapham 1991), in the reaction of catalytic oxidation 
of carbon oxide (Jakubith, Rotermund, Engel, von Oertzen & Ertl 1990), rusting of the steel 
surface in acid with the air (Agladze & Steinbock 2000), in liquid crystal (Frisch, Rica, Coullct 
& Gilli 1994) and laser (Yu, Lu & Harrison 1999) systems. On a larger scale, there are waves of 
infectious diseases travelling through biological populations (Carey, Giles & Mclean 1978, Murray, 
Stanley & Brown 1986), and spiral galaxies (Madore & Freedman 1987, Schulman & Seiden 1986). 
Yet for experimental studies of spiral waves dynamics the BZ reaction medium remains the most 
favourite. 



A common feature of all these phenomena is that they can be mathematically approximated 
by "reaction-diffusion" partial differential equations, 

9 t u = f(u)+DV 2 u, u,fel',Del ,xf ,l>2, (1) 

where u(r, t) is a column- vector of the reagent concentrations, f(u) of the reaction rates, D is 
the matrix of diffusion coefficients, and r e R 2 is the vector of coordinates on the plane. Since 
these equations are essentially nonlinear, their spiral wave solutions in general case are studied 
numerically. Thus, given the complexity of the problem, the current understanding of spiral waves 
is mostly empirical and gives neither possibility for systematic quantitative predictions of the drift , 
nor general understanding on how to control the smooth dynamics of autowave vortices, which is 
important for many practical applications. Effective control of re-entry in excitable cardiac tissue 
will provide a solution to dangerous arrhythmias and fatal fibrillation. 

As a model self-organizing structure, spiral wave demonstrates a remarkable stability, just 
changing its rotational frequency and core location, i.e. drifting, in response to small pertur- 
bations of the medium. As experiments with BZ reaction medium (Agladzc 2000) and computer 
simulations showed spiral waves insensitivity to distant events, it was conjectured (Biktashev 1989) 
that the RFs must decay quickly with distance from the spiral wave core, i.e. spiral waves look 
like essentially non-localized regimes but behave as effectively localized particles (Biktasheva & 
Biktashev 2003). The asymptotical theory of the spiral wave drift, proposed in (Keener 1988, Bik- 
tashev & Holden 1995) and shortly described below, is based on the idea of summation of elemen- 
tary responses of the spiral wave core position and rotation phase to elementary perturbations of 
different modalities and at different times and places. This is mathematically expressed in terms 
of the spiral wave response functions (RFs) equal to zero in the region where the spiral wave is 
insensitive to small perturbations. 

So far, the response functions have been explicitly found with good quantitative accuracy 
only for spiral waves in oscillatory medium described by the Complex Ginzburg-Landau Equa- 
tion (CGLE) (Biktasheva, Elkin & Biktashev 1998, Biktasheva & Biktashev 2001, Biktasheva & 
Biktashev 2003). It was shown that the response functions of vortices in the CGLE medium are 
essentially nonzero only in the vicinity of the core for all sets of model parameters stable spiral 
wave solution exists for (Biktasheva & Biktashev 2001, Biktasheva & Biktashev 2003), which ex- 
plains the localised sensitivity of spiral waves to small perturbations. Most important is the RFs 
ability to make quantitative prediction of spiral wave drift velocity due to small perturbations of 
any nature. 

Thus, a spiral wave organise the medium dividing it into two unequal parts, the core, events 
in which arc translated throughout the medium, and the periphery, obeying the signals from the 
core. It creates a macroscopic wave-particle dualism as an emergent property of the nonlinear 
field, when the regime appears as a non-localized object filling up all available space, but behaves 
as a localized object, only sensitive to perturbations affecting its core. 

Another class of media supporting spiral waves are excitable media. These are even of more 
interest than the oscillatory ones, due to their role in the cardiac, smooth muscle research and 
neuroscience. In order to check localisation properties of the response functions of vortices in 
excitable media, the RFs need to be found explicitly for a particular excitable model. Hamm 
(Hamm 1997) tried to find the response functions for the Barkley model of an excitable system. 
The obtained response functions were effectively localised in the vicinity of the spiral wave core, 
but the accuracy of the solution was not sufficient to allow it to be used for prediction of the 
velocity of the spiral wave drift. 

In this paper we find the response functions for a spiral wave solution in the FitzHugh-Nagumo 
(FHN) model (FitzHugh 1961, Nagumo, Arimoto & Yoshizawa 1962) and show that the RFs are 
effectively localized in the vicinity of the spiral wave core. The model parameters were selected 
to produce an excitable medium with a rigidly rotating spiral wave. The method of computation 
is based on the idea of a moving frame of reference, whose movement is controlled by the spiral 
wave solution found in that frame (Biktashev, Holden & Nikolaev 1996). 

The FitzHugh-Nagumo model is historically the first simplified model of biological excitation. 
It has been studied and used as a classic model for computer simulation of spiral wave dynamics 
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for decades, for it captures the key phenomena of the excitable media while consisting of just two 
partial differential equations, which makes the FHN model easy to study both numerically and 
analytically, 

dtui = e _1 (ui — uf/3 — U2) + D\V 2 u\ 

d t u 2 = e(ui+/3- 7u 2 ) (2) 

where /3, 7, e and D\ are parameters. 

In the FHN model the variable u\ is the fast variable, corresponding to the voltage in biophysi- 
cally realistic models of membrane action potential, and 1*2 does not have any specific physiological 
interpretation, just plays the role of the slow "recovery" variable. The cubic nonlinearity of the 
system results in the simple N-shape nullcline on the phase portrait and explains the key aspects 
of excitability. Existence of spiral wave solutions in the FHN model, their characteristics and 
behaviour depending on the model parameters have been extensively studied by many authors. 
The classic review on the subject is the Winfree article (Winfree 1991). 

2 Asymptotic theory of spiral waves dynamics 

2.1 Initial definitions 

Consider a slightly perturbed "reaction-diffusion" system (1) in two spatial dimensions, 

d t u = f(u) +DV 2 u + eh, u,f,heM f ,DeM (x ',f>2, (3) 

where eh(u, r, t) is a small perturbation, and f el 2 . 

We assume that unperturbed system (1) has solutions in the form of steadily rotating spiral 
waves, 

u = V(f,t)=V(p(r),'d(f f )+ujt). (4) 

Here 

6 = tf(r) + cut 

is a polar angle in the "corotating" frame of reference, which is rigidly rotating with the angular 
velocity cj, while p(f) and i?(f) are the polar coordinates in the original (laboratory) frame of 
reference. 

The unperturbed reaction-diffusion system (1), or (3) with eh = 0, has an obvious but impor- 
tant symmetry: it is invariant with respect to the Euclidean group of motions of the plane {r}. 
Since solution (4) at any fixed t is not invariant against this group, the group "multiplies" this 
solution. That is, 

U'(f, t) = U(p(f - &),<&(?- R) + Q), (5) 

where = ujt — <!>, is another solution for any constant displacement vector R = (X, Yy and 
initial rotation phase $. 

Thus, if the unperturbed system has one spiral wave solution, then it has a whole three- 
dimensional manifold of such solutions, that are relatively stable with respect to the shift along 
the manifold. 

2.2 Finite-dimensional analogy 

The asymptotic theory of drift of spiral waves (Biktashcv & Holden 1995) was proposed based on 
the analogy with finite-dimension problem of perturbation of an invariant manifold (sec fig. 1). 
If a vector field f(u) in an n-dimcnsional phase space has an invariant m-dimensional manifold 
U(a), m < n, stable as a whole, then small perturbation of this vector field will, under certain 
conditions, preserve the invariant manifold, just slightly displacing it, U >—>■ U'. Another effect 
of the perturbation is that the vector field on the shifted manifold A' (a) will be slightly different 
from the original one, A (a). In practice, the existence of the original invariant manifold U(a) 
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Figure 1: Perturbation of an invariant manifold. Vector field f (u) in phase space with coordinates 
u has an invariant manifold U with coordinates a, and vector field A on the manifold. Perturbed 
vector field f'(u) has a slightly different invariant manifold, U', and a slightly different vector field 
A' on it. Original objects are shown by solid lines, and perturbed objects by dashed lines. 



could be due to a symmetry group. In that case, the flow on that manifold could be in some sense 
degenerate, and then the perturbation will remove this degeneracy. 

To compare the two vector fields, on the original manifold and on the perturbed, we need to 
relate their coordinate systems {a}. A natural way is to require that the vector connecting two 
corresponding points U(a) and U'(a), would not have a component along the manifold, i.e. along 
any of the tangent vectors Vj(a) = dXJ/daj. In other words, it should be orthogonal, 

(W J (a),U , (a)-U(a)} = 0, j = l...m, (6) 

to the projectors Wj(a) onto the tangent vectors Vj(a): 

(W i (o),V fc (o)>=^, fc . (7) 

These projectors are eigenvectors of the adjoint linearized matrix (di / du) T (a) . The two effects of 
the perturbation are produced by its two components, along and across the manifold, as determined 
by the projectors W 3 -. 

Thus, if the manifold comprises only non-moving points, the tangent component will determine 
the slow drift along the manifold. 

If this finite-dimensional scheme can be applied to spiral waves, the role of the vector field is 
played by the reaction-diffusion system, so the phase space is a functional space. The invariant 
manifold is the three-dimensional manifold of spiral waves and is due to a symmetry group, the 
Euclidean group of the plane. The coordinates on the manifold are ReK 2 , the centre of rotation 
of the spiral wave, and 0, its rotation angle. The flow on the manifold is degenerate, as it consists 
of relatively stable periodic orbits, which correspond to steady rotation of spiral waves around 
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fixed centres: 

9 = u>t - <f>, <f> = const; R = const. (8) 

The perturbation removes this degeneracy, and we observe the drift of the spirals. By analogy 
with the finite-dimensional case, we expect that the flow on the manifold of spiral waves will be 
described by 

d t Q = LJ + eH (R,@), d t R = sHx(R,e), (9) 

where Hq and Hi are "projections" of the perturbation onto the tangent space of the manifold 
U(a). The right-hand sides of (9) depend on the phase 9. 9n the time scale £ _1 this phase 
oscillates fast; averaging over these oscillations gives motion equations of the spiral waves, 

d t e = uj + sH^(R) + 0(e 2 ), d t R= ei?i (R) + O (e 2 ) . (10) 
2.3 Response functions 

Thus, the finite-dimensional analogy suggests that the dynamics of spiral waves (perhaps like that 
of many other dissipative structures) is described by "Aristotelean" mechanics, when the velocity 
of motion is proportional to the applied perturbation. The right-hand sides in the equations, 
the "forces" , are projections of the perturbation onto the corresponding tangent space of the 
invariant manifold U(a). This tangent space is a linear space, the span of the Goldstonc modes, 
corresponding to the translations along the symmetry group, at R = and 9 = 0, 



V = -^- 1 a t U(r,t)| t=o = -a e U(p(r),0(r)), 
V±i = {d x Tid v )WMt=o = ~ {®x T idy) U(r, t) = - l -e^ e (d^ip^de) fj(p(f),6(r)) 

(11) 



Here mode Vo corresponds to the shift in time (or to what is the same, rotation in space), and 
Vi corresponds to the shift in space. Tildes in (11) designate the corotating frame of reference, 
so x, y are Cartesian coordinates there and U is the unperturbed spiral wave solution, which is 
stationary in that frame of reference. We omit the tildes henceforth for brevity. 
The Goldstone modes are critical cigenfunctions 

CV n = iunV n , n = 0, ±1 (12) 

of the linearized operator C: 



Oil 



(13) 



Here again the tilde at C reminds that the linearized operator is considered in the co-rotating 
frame of reference where it does not depend on time. The additive tl —u)d$" appears here due to 
rotation with respect to the original system of coordinates. 

Thus, for each particular point at the manifold, the projection operators map the functional 
space of the perturbations into the three-dimensional tangent space, and are thus just three func- 
tionals. Since all points of our manifold are equivalent to each other up to a Euclidean transfor- 
mation of the plane, it is enough to know the projection functionals at one point, i.e. just for one 
location of the spiral wave. This symmetry consideration shows that if the functionals H n are 
written as integrals, they should have the form: 

H^) = e m * { — [ [ d 2 fe' m ^ T (W n (p(r- R),^- R)+lut -<i>) ,h) , (14) 



2tt 

t — 7T / LU 
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where 



h = h(U(r,T),r,T), 
R = R(t), 
$ = §(*), 

ST=(li) a +i(li)v (15) 
The kernels W„ of the integrals (14) are eigenfunctions 

£ + W n = -iu;nW n , n = 0, ±1. (16) 

of the adjoint linearized operator considered in the co-rotating frame of reference: 

(17) 

u=U(r) 

As in the finite dimensional example, we assume here that W are normalized in such a way that 
(7) is satisfied. 

The functions W„ are called the response functions (RFs) of the spiral wave. Folloring the 
analogy with the Goldstone modes, Wo defines the shift in time (or the turning in space), and 
Wi defines the shift in space. So Wo is called the temporal or rotation RF, and Wi is called the 
spatial or shift RF. 



C+ = DV 2 + code + 




3 Mathematical formulation of the problem for the FitzHugh- 
Nagumo model 

The problems (1), (4) and (16)-(17) for the FitzHugh-Nagumo system (2) take the form 
e" 1 (Ut - Uf/i - U 2 ) + (.DiV 2 - ud e )U x = 0, 

e(JJx+P-lU i )-wdt)U 2 = Q, (18) 

(e-^l - Ul) + w(in + do) + ^iV 2 ) W{ 1 + eW% = 0, 

-e _1 W7 + (-ey + u(in + d g )) W' 2 l = 0, n = 0, 1 (19) 

for the unperturbed solution Uj, its angular velocity u) and the RFs W™, n = 0,1, j = 1,2, 
€ M, W12 6 C • System (18)— (19) should be supplied with normalisation and boundary 
conditions, and discretized. For discretisation we used rectangular grids in Cartesian coordinates. 

3.1 Spiral wave problem 

Since the Response Functions are the solution of the adjoint linearized problem in the system of 
reference that is rotating with the angular velocity of the spiral itself, we need first to find the 
spiral wave solution in this system of reference, i.e. we need to find the spiral wave solution together 
with its rotation angular velocity. So we have a nonlinear eigenvalue together with a boundary 
value problem. 

We solved this problem numerically on a square domain [x, y) G S = [-L/2, L/2] x [-L/2, L/2], 
for different L from 25 to 50. First, the spiral wave was initiated by solving a Cauchy problem 
for (2) for initial conditions u(x,y,0) = 0.7sign(x), v(x,y,0) — 0.6sign(y). When a stationary 
rotating spiral wave was established, typically within time interval t £ [0, T], T ~ 40, the resulting 
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distribution u(x,y,T), v(x,y,T) was used as an initial condition for the following system: 
dtux = e^ 1 (u 1 ~ul/3-u 2 )+D 1 V 2 u 1 + ^ CjdjUx, 

j=x,y,8 

d t u 2 = e (iti + [3 - 7U2) + CjdjU 2 , 

j=x,y,6 

Cj = -q 3 C } ■- ^ J {dtuxdjui + d t u 2 djU 2 )dxdy, j = x,y,9, (20) 
v 

where de = xd y - yd x , V = {(x, y) : x 2 + y 2 < (L/2) 2 }, A = ir(L/2) 2 , q e = 0, and positive 
coefficients p x ,y,e and q XiV have been selected by experimentation. Informally, the idea of this 
system, adopted with appropriate modification from (Biktashev et al. 1996), is that the first two 
equations are system (2) in a frame of reference moving with speeds —C XiV $ in the x, y and 8 
directions, and the integrals in the evolution equations for C x ^ v> g are "detectors of movement" in 
those directions. So the movement of the frame of reference is adjusted in such a way so as to make 
the solution in this frame of reference is stationary. On the formal level, it is straightforward to see 
that if solution of (20) converges to a stationary state, then ui l2 will satisfy (18) with u> = —Cg, 
neglecting the boundary conditions. 

We considered the problem (20) for the following set of parameters: D± = 1.0, e = 0.30, 
j3 = 0.75, 7 = 0.50, which correspond to a rigidly rotating spiral wave solution (Winfree 1991). 
Calculations were performed using the explicit Euler method in time, central differences in space, 
with fixed time step from At = 3 • 10~ 3 down to At = 5 • 10~ 4 and space step Aa; = 0.5, with 
Neuman boundary conditions on a rectangular grid in Cartesian coordinates: d x ui(±L/2,y) = 
d y ui(x, ±L/2) = 0. The frame of reference adjustment parameters were chosen q x ^ y = 1, p x ^ y = 7 
and pe = 5. 

The result is a stable, stationary spiral (demonstrated on fig. 2) in the system of reference, 
which rotates with the angular velocity u> « 0.32 clockwise. Having this angular velocity and the 
spiral wave solution, it is possible to find the corresponding response functions. 

3.2 The response functions problem 

The response functions W were calculated simultaneously with finding the spiral wave solution 
for (20), by solving the adjoint linearized problems 

dtWi = e _1 (l — uf)wi + ew 2 + DV 2 wi — CjdjWi, 

j=x,y,9 

dtw 2 = -e~ 1 w 1 - e~fw 2 - ^2 C 3 djW 2 . (21) 

j=x,y,9 

As (ui, u 2 , Ci, C 2 , C3), solution of (20), converges to (Ui, U 2 , 0, 0, — u>), solution of (19), then w = 
(u>i, w 2 ) T , a typical solution of (21), is expected to converge to 

w(x, y, t) a cdW (i, y) + Re (dWi(x, y)e~ iMt ) , (22) 

where cq G R and c\ E C, ideally, are constants depending on initial conditions, but in real calcu- 
lations can slowly change in time due to numerical approximation. To obtain Wo,i, we calculate 
three solutions w'™), m = 1,2,3, to (21) with three linearly independent initial conditions. As- 
suming that, after a sufficiently long time, these satisfy (22), we can take their appropriate linear 
combinations to satisfy (7). So for a given triplet w' m \ m = 1,2,3, we are looking for constants 
Pj,m, j = $ 5 y, m = 1, 2, 3, so that 

m—1 
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Figure 2: Spiral wave solution and the response functions. Parameters: L = 50, Ax = 0.5, 
At = 5 • 10 -4 , Neuman boundary conditions on dS for U and Dirichlet boundary conditions on 
dV for W. 
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would satisfy biorthogonality condition (7) with respect to the set Vft = <9fcU, k = 6,x,y. This 
requirement implies that 

3 

£p,- TO (wW U fe ) = £ jVbl 

m— 1 

i.e. P = (Pj,m) is an inverse to the matrix B = (B„ u k), where B m ^ = (w^" l ',Ufe) are obtained 
by integrating the given solutions of (21) with the derivatives of the solution of the nonlinear 
problem. Having thus found W 1jS) b after integrating (21) for a time interval of certain r, we used 
these W li9 .8 as three linearly independent initial conditions for (21) for a further time interval of 
length t. 

Equation (21) was integrated on the disk (x,y) G T> C S with Dirichlet boundary conditions 
wi(x,y) = W2(x 1 y) = 0, [x, y) G dT> , which were implemented by setting w\ t 2{x,y) = for 
(x, y) G S \ V. We tried r = 1 and r = 10 with identical results; the solution obtained by this 
procedure converged within time scale of t ~ 20. As Vft = <9fcU are related to the Goldstone 
modes Vo,±i by (11), this gives relationship between functions ^ x ,y,e found in this way, with 
RFs, in the form W = - W e and Wj = -{W x - iW y ). 

Figure 2 shows the components of the Response Functions obtained in this way. The most 
important result is that all components are localized in a close vicinity of the tip of the spiral. It 
can also be observed that the amplitude of the second components of the RFs is higher than that 
of the first, so controlling movement of the spiral wave via the second component, the inhibitor, 
if it is pratically feasible, could be more efficient than via the first component the activator. 

The accuracy of the solution has been checked by varying the discretization steps At and Ax, 
the domain size L, boundary conditions (Neuman vs Dirichlet, dT> vs dS). Based on such checks, 
we estimate that the accuracy of the solutions is within a few percent. Further improvement of ac- 
curacy requires decrease of Ax while keeping L same or increasing, and due to stability limitations 
of the fully explicit Euler time stepping, such improvement is extremely costly if done within the 
same numeric scheme, and requires a more advanced scheme, e.g. fully implicit/pscudospectral in 
the 9 direction. 

4 Conclusion 

The response functions are very important characteristics of the spiral wave, for they define the 
phenomenology of the spiral behaviour. Experiments and computer simulations, demonstrating 
the spirals insensitivity to distant events, implied that the vortices Response Functions must 
decay quickly with distance from the core. Such decay will guarantee the convergence of the 
integrals (14) even for non-localized perturbations, for example caused by the simultaneous change 
of the properties of the whole medium. But the mathematical peculiarity of the idea presuming 
qualitatively different behaviour of eigenfunctions of the linear operator and its adjoint one resulted 
in a natural scepticism. 

Recently shown for the CGLE oscillatory medium localisation of the response functions in the 
vicinity of the vortices core left open the question on localisation properties of the spiral waves 
response functions in an excitable medium. The results obtained in this paper confirm the exis- 
tence of effectively localized response functions of a spiral wave solution in the classical excitable 
FitzHugh-Nagumo model, at least for the particular set of the model parameters corresponding 
to a rigidly rotating spiral wave. 
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